Laplace distribution: laplace (double exponential)

Laplace distribution: laplace (double exponential)#

The Laplace distribution is a symmetric, continuous distribution with a sharp peak at its location and exponentially decaying tails. Compared to a Gaussian, it puts more mass near the center and more mass in the tails, which makes it a common choice for robust modeling.

In SciPy it appears as scipy.stats.laplace.

Learning goals#

  • understand what the Laplace distribution models and when it is useful

  • write down the PDF/CDF and connect them to sampling and likelihood

  • compute key moments (mean/variance/skewness/kurtosis) and entropy

  • derive the closed-form MLE (median + mean absolute deviation)

  • implement NumPy-only sampling and validate everything by Monte Carlo

Notebook roadmap#

  1. Title & Classification

  2. Intuition & Motivation

  3. Formal Definition

  4. Moments & Properties

  5. Parameter Interpretation

  6. Derivations

  7. Sampling & Simulation (NumPy-only)

  8. Visualization (PDF/CDF + Monte Carlo)

  9. SciPy Integration

  10. Statistical Use Cases

  11. Pitfalls

  12. Summary

import numpy as np
import plotly
import plotly.express as px
import plotly.graph_objects as go
import os
import plotly.io as pio
from plotly.subplots import make_subplots

import scipy
from scipy.stats import chi2, laplace, norm

pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")

np.set_printoptions(precision=4, suppress=True)

SEED = 7
rng = np.random.default_rng(SEED)

print("NumPy ", np.__version__)
print("SciPy ", scipy.__version__)
print("Plotly", plotly.__version__)

NumPy  1.26.2
SciPy  1.15.0
Plotly 6.5.2

1) Title & Classification#

  • Name: laplace

  • Type: continuous distribution

  • Support: x ∈ (-∞, ∞)

  • Parameter space: location μ ∈ ℝ and scale b > 0

We write:

\[X \sim \mathrm{Laplace}(\mu, b).\]

The standard Laplace is (\mathrm{Laplace}(0, 1)).

2) Intuition & Motivation#

What it models#

The Laplace distribution is a natural model for real-valued errors that are:

  • centered around a typical value (μ)

  • more sharply peaked than a Gaussian

  • heavy-tailed relative to a Gaussian, but still with exponential tails

A clean way to remember it:

  • Normal: log-density is quadratic in ((x-\mu))

  • Laplace: log-density is linear in (|x-\mu|)

So Laplace errors correspond to an L1 loss.

Typical real-world use cases#

  • Robust regression / signal processing: modeling residuals with Laplace leads to least-absolute-deviation fitting.

  • Sparsity-inducing priors: the Laplace prior underpins the Bayesian view of Lasso (L1 regularization).

  • Differential privacy (Laplace mechanism): adding Laplace noise to numeric queries for privacy guarantees.

  • Noise with occasional big jumps: exponential tails can be a better match than Gaussian in some domains.

Relations to other distributions#

  • If (E_1, E_2) are i.i.d. exponential with mean (b), then (E_1 - E_2 \sim \mathrm{Laplace}(0, b)).

  • If (S\in{-1,+1}) is a fair sign and (E\sim\mathrm{Exp}(\text{mean}=b)), then (S,E\sim\mathrm{Laplace}(0, b)).

  • If (X\sim\mathrm{Laplace}(\mu, b)), then (|X-\mu|\sim\mathrm{Exp}(\text{mean}=b)).

  • The Laplace has heavier tails than a Gaussian but all moments exist; its MGF exists only for (|t|<1/b).

3) Formal Definition#

PDF#

For (X\sim\mathrm{Laplace}(\mu, b)):

\[ f(x; \mu, b) = \frac{1}{2b}\,\exp\left(-\frac{|x-\mu|}{b}\right), \qquad b>0. \]

CDF#

The CDF has a simple piecewise form:

\[\begin{split} F(x; \mu, b) = \begin{cases} \tfrac{1}{2}\exp\left(\frac{x-\mu}{b}\right), & x<\mu\\ 1 - \tfrac{1}{2}\exp\left(-\frac{x-\mu}{b}\right), & x\ge \mu. \end{cases} \end{split}\]

Quantile function (inverse CDF)#

For (p\in(0,1)):

\[\begin{split} F^{-1}(p) = \begin{cases} \mu + b\,\log(2p), & 0<p<\tfrac{1}{2}\\ \mu - b\,\log\big(2(1-p)\big), & \tfrac{1}{2}\le p<1. \end{cases} \end{split}\]
def _check_scale(b: float) -> float:
    b = float(b)
    if not np.isfinite(b) or b <= 0:
        raise ValueError("`b` (scale) must be a positive, finite number.")
    return b


def laplace_pdf(x, mu: float = 0.0, b: float = 1.0) -> np.ndarray:
    b = _check_scale(b)
    x = np.asarray(x, dtype=float)
    return np.exp(-np.abs(x - mu) / b) / (2.0 * b)


def laplace_logpdf(x, mu: float = 0.0, b: float = 1.0) -> np.ndarray:
    b = _check_scale(b)
    x = np.asarray(x, dtype=float)
    return -np.log(2.0 * b) - np.abs(x - mu) / b


def laplace_cdf(x, mu: float = 0.0, b: float = 1.0) -> np.ndarray:
    b = _check_scale(b)
    x = np.asarray(x, dtype=float)
    z = (x - mu) / b
    return np.where(x < mu, 0.5 * np.exp(z), 1.0 - 0.5 * np.exp(-z))


def laplace_ppf(p, mu: float = 0.0, b: float = 1.0, eps: float = 1e-12) -> np.ndarray:
    b = _check_scale(b)
    p = np.asarray(p, dtype=float)
    if np.any((p <= 0) | (p >= 1)):
        raise ValueError("p must be in (0, 1)")
    p = np.clip(p, eps, 1.0 - eps)
    left = mu + b * np.log(2.0 * p)
    right = mu - b * (np.log(2.0) + np.log1p(-p))  # log(2(1-p)) stably
    return np.where(p < 0.5, left, right)


# Quick cross-check vs SciPy
mu, b = 0.4, 1.7
x = np.linspace(-3, 3, 9)
p = np.array([0.1, 0.5, 0.9])

print("max |pdf - scipy|:", float(np.max(np.abs(laplace_pdf(x, mu, b) - laplace.pdf(x, loc=mu, scale=b)))))
print("max |cdf - scipy|:", float(np.max(np.abs(laplace_cdf(x, mu, b) - laplace.cdf(x, loc=mu, scale=b)))))
print("max |ppf - scipy|:", float(np.max(np.abs(laplace_ppf(p, mu, b) - laplace.ppf(p, loc=mu, scale=b)))))

max |pdf - scipy|: 0.0
max |cdf - scipy|: 0.0
max |ppf - scipy|: 0.0

4) Moments & Properties#

For (X\sim\mathrm{Laplace}(\mu, b)):

  • Mean: (\mathbb{E}[X]=\mu)

  • Variance: (\mathrm{Var}(X)=2b^2) (so (\mathrm{sd}(X)=\sqrt{2},b))

  • Skewness: 0 (symmetry)

  • Kurtosis: 6 (Pearson), so excess kurtosis = 3

MGF and characteristic function#

  • MGF exists only for (|t|<1/b):

\[ M_X(t) = \mathbb{E}[e^{tX}] = \frac{e^{\mu t}}{1-b^2 t^2}, \qquad |t|<1/b. \]

  • Characteristic function (exists for all real (t)):

\[ \varphi_X(t) = \mathbb{E}[e^{itX}] = \frac{e^{i\mu t}}{1+b^2 t^2}. \]

Entropy#

The differential entropy is:

\[ H(X) = 1 + \log(2b). \]
def laplace_moments(mu: float = 0.0, b: float = 1.0):
    b = _check_scale(b)
    mean = float(mu)
    var = float(2.0 * b * b)
    skew = 0.0
    kurt_excess = 3.0
    entropy = float(1.0 + np.log(2.0 * b))
    return mean, var, skew, kurt_excess, entropy


def sample_moments(x: np.ndarray):
    x = np.asarray(x, dtype=float)
    m = float(x.mean())
    c = x - m
    v = float(np.mean(c**2))
    skew = float(np.mean(c**3) / (v ** 1.5))
    kurt_excess = float(np.mean(c**4) / (v**2) - 3.0)
    return m, v, skew, kurt_excess


mu, b = 1.5, 0.8

mean_f, var_f, skew_f, kurt_f, ent_f = laplace_moments(mu, b)
mean_s, var_s, skew_s, kurt_s = laplace.stats(loc=mu, scale=b, moments="mvsk")
ent_s = laplace.entropy(loc=mu, scale=b)

print("theory  (mean, var, skew, kurt_excess, entropy):", (mean_f, var_f, skew_f, kurt_f, ent_f))
print("scipy   (mean, var, skew, kurt_excess, entropy):", (float(mean_s), float(var_s), float(skew_s), float(kurt_s), float(ent_s)))

x = laplace.rvs(loc=mu, scale=b, size=300_000, random_state=rng)
mean_mc, var_mc, skew_mc, kurt_mc = sample_moments(x)
print("monte   (mean, var, skew, kurt_excess):         ", (mean_mc, var_mc, skew_mc, kurt_mc))

theory  (mean, var, skew, kurt_excess, entropy): (1.5, 1.2800000000000002, 0.0, 3.0, 1.4700036292457357)
scipy   (mean, var, skew, kurt_excess, entropy): (1.5, 1.2800000000000002, 0.0, 3.0, 1.4700036292457357)
monte   (mean, var, skew, kurt_excess):          (1.501420070672676, 1.2685335280033816, -0.0003665014939261965, 2.9428857779033946)

def laplace_cf(t, mu: float = 0.0, b: float = 1.0) -> np.ndarray:
    b = _check_scale(b)
    t = np.asarray(t, dtype=float)
    return np.exp(1j * mu * t) / (1.0 + (b * t) ** 2)


mu, b = 0.0, 1.0
t = np.linspace(-15, 15, 2000)
phi = laplace_cf(t, mu=mu, b=b)

fig = make_subplots(rows=1, cols=2, subplot_titles=("Re φ(t)", "Im φ(t)"))
fig.add_trace(go.Scatter(x=t, y=np.real(phi), mode="lines"), row=1, col=1)
fig.add_trace(go.Scatter(x=t, y=np.imag(phi), mode="lines"), row=1, col=2)
fig.update_xaxes(title_text="t", row=1, col=1)
fig.update_xaxes(title_text="t", row=1, col=2)
fig.update_layout(width=950, height=350, showlegend=False, title="Characteristic function of Laplace(0,1)")
fig.show()

5) Parameter Interpretation#

  • μ (location) shifts the distribution left/right. For Laplace it is simultaneously the mean, median, and mode.

  • b (scale) controls both the peak height and tail thickness:

    • (f(\mu)=1/(2b)) so larger b means a lower peak.

    • (\mathrm{sd}(X)=\sqrt{2},b) so b is proportional to standard deviation.

Below we visualize how ((\mu,b)) changes the PDF and CDF.

x = np.linspace(-10, 10, 3000)
params = [
    (0.0, 0.5),
    (0.0, 1.0),
    (0.0, 2.0),
    (2.0, 1.0),
]

fig = make_subplots(rows=1, cols=2, subplot_titles=("PDF", "CDF"))
for mu, b in params:
    fig.add_trace(
        go.Scatter(x=x, y=laplace_pdf(x, mu=mu, b=b), mode="lines", name=f"μ={mu}, b={b}"),
        row=1,
        col=1,
    )
    fig.add_trace(
        go.Scatter(x=x, y=laplace_cdf(x, mu=mu, b=b), mode="lines", name=f"μ={mu}, b={b}"),
        row=1,
        col=2,
    )

fig.update_xaxes(title_text="x", row=1, col=1)
fig.update_yaxes(title_text="f(x)", row=1, col=1)
fig.update_xaxes(title_text="x", row=1, col=2)
fig.update_yaxes(title_text="F(x)", row=1, col=2)
fig.update_layout(width=1000, height=420)
fig.show()

6) Derivations#

Expectation#

Let (Y=X-\mu). Then (Y) has symmetric density (f_Y(y) = (1/(2b))\exp(-|y|/b)).

Because (y f_Y(y)) is an odd function,

\[ \mathbb{E}[Y] = \int_{-\infty}^{\infty} y f_Y(y)\,dy = 0 \quad\Rightarrow\quad \mathbb{E}[X] = \mu. \]

Variance#

Using symmetry again,

\[ \mathbb{E}[(X-\mu)^2] = \int_{-\infty}^{\infty} y^2\frac{1}{2b}e^{-|y|/b}dy = \frac{1}{b}\int_{0}^{\infty} y^2 e^{-y/b}dy. \]

The integral is a Gamma-type integral: (\int_0^{\infty} y^2 e^{-y/b}dy = 2!,b^3 = 2b^3), hence

\[\mathrm{Var}(X)=\mathbb{E}[(X-\mu)^2]=2b^2.\]

Likelihood (and the MLE)#

For i.i.d. data (x_1,\dots,x_n) from (\mathrm{Laplace}(\mu,b)), the likelihood is

\[ L(\mu,b) = \prod_{i=1}^n \frac{1}{2b}\exp\left(-\frac{|x_i-\mu|}{b}\right) = (2b)^{-n}\exp\left(-\frac{1}{b}\sum_{i=1}^n |x_i-\mu|\right). \]

So the log-likelihood is

\[ \ell(\mu,b) = -n\log(2b) - \frac{1}{b}\sum_{i=1}^n |x_i-\mu|. \]

Key consequence: for fixed (b), maximizing (\ell) over (\mu) is equivalent to minimizing (\sum_i |x_i-\mu|), which is minimized by any median of the sample.

With (\hat\mu) chosen as a sample median, differentiating (\ell(\hat\mu,b)) with respect to (b) gives the (closed-form) MLE:

\[\hat b = \frac{1}{n}\sum_{i=1}^n |x_i-\hat\mu|.\]
def laplace_loglik(data: np.ndarray, mu: float, b: float) -> float:
    return float(np.sum(laplace_logpdf(data, mu=mu, b=b)))


def laplace_mle_closed_form(data: np.ndarray, b_floor: float = 1e-12) -> tuple[float, float]:
    """Closed-form MLE for Laplace(μ,b): μ=median, b=mean absolute deviation from μ.

    Notes
    -----
    - If all observations are identical, the likelihood increases as b -> 0+.
      We return a small positive floor to stay inside the parameter space.
    """
    x = np.asarray(data, dtype=float)
    if x.ndim != 1:
        raise ValueError("data must be 1D")
    mu_hat = float(np.median(x))
    b_hat = float(np.mean(np.abs(x - mu_hat)))
    b_hat = float(max(b_hat, b_floor))
    return mu_hat, b_hat


mu_true, b_true = 1.0, 0.7
data = laplace.rvs(loc=mu_true, scale=b_true, size=2000, random_state=rng)

mu_hat, b_hat = laplace_mle_closed_form(data)
mu_hat_sp, b_hat_sp = laplace.fit(data)

print("true (μ, b)         =", (mu_true, b_true))
print("closed-form MLE     =", (mu_hat, b_hat))
print("scipy laplace.fit   =", (float(mu_hat_sp), float(b_hat_sp)))
print("loglik at MLE       =", laplace_loglik(data, mu_hat, b_hat))

true (μ, b)         = (1.0, 0.7)
closed-form MLE     = (1.000653645441159, 0.716006489989578)
scipy laplace.fit   = (1.000653645441159, 0.716006489989578)
loglik at MLE       = -2718.1622654572566

7) Sampling & Simulation (NumPy-only)#

Two convenient sampling views:

  1. Inverse CDF: sample (U\sim\mathrm{Uniform}(0,1)) and return (F^{-1}(U)).

  2. Difference of exponentials: if (E_1,E_2\overset{iid}{\sim}\mathrm{Exp}(\text{mean}=b)), then (\mu + (E_1-E_2)\sim\mathrm{Laplace}(\mu,b)).

We’ll implement both with NumPy and sanity-check them.

def laplace_rvs_inverse(
    rng: np.random.Generator,
    size: int,
    mu: float = 0.0,
    b: float = 1.0,
    eps: float = 1e-12,
) -> np.ndarray:
    b = _check_scale(b)
    u = rng.random(size)
    u = np.clip(u, eps, 1.0 - eps)
    left = mu + b * np.log(2.0 * u)
    right = mu - b * (np.log(2.0) + np.log1p(-u))
    return np.where(u < 0.5, left, right)


def laplace_rvs_exp_difference(
    rng: np.random.Generator,
    size: int,
    mu: float = 0.0,
    b: float = 1.0,
) -> np.ndarray:
    b = _check_scale(b)
    e1 = rng.exponential(scale=b, size=size)
    e2 = rng.exponential(scale=b, size=size)
    return mu + (e1 - e2)


mu, b = -0.5, 1.2
n = 200_000

x_inv = laplace_rvs_inverse(rng, size=n, mu=mu, b=b)
x_diff = laplace_rvs_exp_difference(rng, size=n, mu=mu, b=b)

q = [0.05, 0.25, 0.5, 0.75, 0.95]
print("theory quantiles:", np.round(laplace.ppf(q, loc=mu, scale=b), 4))
print("inv   quantiles:", np.round(np.quantile(x_inv, q), 4))
print("diff  quantiles:", np.round(np.quantile(x_diff, q), 4))

theory quantiles: [-3.2631 -1.3318 -0.5     0.3318  2.2631]
inv   quantiles: [-3.2641 -1.3348 -0.5003  0.3281  2.2367]
diff  quantiles: [-3.2562 -1.3277 -0.4983  0.3328  2.2574]

8) Visualization (PDF/CDF + Monte Carlo)#

We’ll plot the theoretical PDF/CDF and compare them to Monte Carlo samples.

mu, b = 0.0, 1.0
x = np.linspace(-8, 8, 4000)

fig = make_subplots(rows=1, cols=2, subplot_titles=("PDF", "CDF"))
fig.add_trace(go.Scatter(x=x, y=laplace_pdf(x, mu=mu, b=b), mode="lines", name="pdf"), row=1, col=1)
fig.add_trace(go.Scatter(x=x, y=laplace_cdf(x, mu=mu, b=b), mode="lines", name="cdf"), row=1, col=2)
fig.update_xaxes(title_text="x", row=1, col=1)
fig.update_yaxes(title_text="f(x)", row=1, col=1)
fig.update_xaxes(title_text="x", row=1, col=2)
fig.update_yaxes(title_text="F(x)", row=1, col=2)
fig.update_layout(width=950, height=380, showlegend=False, title="Laplace(0,1): PDF and CDF")
fig.show()


# Monte Carlo histogram vs theory
n = 120_000
samples = laplace_rvs_inverse(rng, size=n, mu=mu, b=b)

fig = go.Figure()
fig.add_trace(
    go.Histogram(
        x=samples,
        histnorm="probability density",
        nbinsx=120,
        name="Monte Carlo",
        opacity=0.6,
    )
)
fig.add_trace(go.Scatter(x=x, y=laplace_pdf(x, mu=mu, b=b), mode="lines", name="Theory"))
fig.update_layout(
    title="Laplace(0,1): histogram vs PDF",
    xaxis_title="x",
    yaxis_title="density",
    width=900,
    height=420,
    barmode="overlay",
)
fig.show()


# Empirical CDF vs theory
xs = np.sort(samples)
ecdf = np.arange(1, xs.size + 1) / xs.size

fig = go.Figure()
fig.add_trace(go.Scatter(x=xs, y=ecdf, mode="lines", name="Empirical"))
fig.add_trace(go.Scatter(x=x, y=laplace_cdf(x, mu=mu, b=b), mode="lines", name="Theory"))
fig.update_layout(title="Laplace(0,1): empirical CDF", xaxis_title="x", yaxis_title="F(x)", width=900, height=420)
fig.show()


# Q-Q plot (sample quantiles vs theoretical quantiles)
p = np.linspace(0.01, 0.99, 200)
q_theory = laplace_ppf(p, mu=mu, b=b)
q_emp = np.quantile(samples, p)

fig = go.Figure()
fig.add_trace(go.Scatter(x=q_theory, y=q_emp, mode="markers", name="quantiles"))
fig.add_trace(go.Scatter(x=q_theory, y=q_theory, mode="lines", name="y=x"))
fig.update_layout(title="Laplace(0,1): Q-Q plot", xaxis_title="theoretical quantile", yaxis_title="sample quantile", width=700, height=520)
fig.show()

9) SciPy Integration#

SciPy’s laplace uses the same parameterization:

  • loc = μ

  • scale = b

The most common methods:

  • laplace.pdf(x, loc, scale) / laplace.logpdf(...)

  • laplace.cdf(x, loc, scale) / laplace.ppf(p, loc, scale)

  • laplace.rvs(loc, scale, size, random_state=...)

  • laplace.fit(data) (MLE)

mu, b = 0.7, 1.2
x = np.array([-2.0, 0.0, 1.0, 3.0])

print("pdf:", laplace.pdf(x, loc=mu, scale=b))
print("cdf:", laplace.cdf(x, loc=mu, scale=b))
print("rvs:", laplace.rvs(loc=mu, scale=b, size=5, random_state=rng))

data = laplace.rvs(loc=mu, scale=b, size=3000, random_state=rng)
mu_hat_sp, b_hat_sp = laplace.fit(data)
mu_hat_cf, b_hat_cf = laplace_mle_closed_form(data)

print("fit (scipy)      :", (float(mu_hat_sp), float(b_hat_sp)))
print("fit (closed-form):", (mu_hat_cf, b_hat_cf))

pdf: [0.0439 0.2325 0.3245 0.0613]
cdf: [0.0527 0.279  0.6106 0.9265]
rvs: [1.3046 0.8555 4.5273 1.1789 0.7094]
fit (scipy)      : (0.7282358708034962, 1.1857197356240277)
fit (closed-form): (0.7282358708034962, 1.1857197356240277)

10) Statistical Use Cases#

Hypothesis testing#

The Laplace likelihood leads to simple likelihood-ratio tests (LRTs). We’ll demonstrate an LRT for the location parameter:

  • (H_0: \mu = \mu_0)

  • (H_1: \mu) free

We keep (b) unknown in both models (estimated by MLE). Asymptotically, the LRT statistic is (\chi^2) with 1 degree of freedom.

Bayesian modeling#

A Laplace prior (p(\theta) \propto \exp(-|\theta|/b)) induces sparsity/shrinkage. Under a Gaussian likelihood, the MAP estimator becomes soft-thresholding (the 1D analog of Lasso).

Generative modeling#

Laplace noise is used as a primitive generator for privacy-preserving releases (Laplace mechanism), and also as a building block in mixture models and robust generative pipelines.

def laplace_lrt_mu(data: np.ndarray, mu0: float) -> tuple[float, float]:
    """Likelihood-ratio test for H0: μ=mu0 vs H1: μ free (b unknown in both)."""
    x = np.asarray(data, dtype=float)
    mu_hat, b_hat = laplace_mle_closed_form(x)

    b0_hat = float(np.mean(np.abs(x - mu0)))
    b0_hat = float(max(b0_hat, 1e-12))

    ll_hat = laplace_loglik(x, mu=mu_hat, b=b_hat)
    ll_0 = laplace_loglik(x, mu=mu0, b=b0_hat)

    LR = 2.0 * (ll_hat - ll_0)
    p_value = float(chi2.sf(LR, df=1))
    return float(LR), p_value


mu0 = 0.0
b_true = 1.0
n = 120

# Under H0
data0 = laplace_rvs_inverse(rng, size=n, mu=mu0, b=b_true)
LR0, p0 = laplace_lrt_mu(data0, mu0=mu0)
print(f"H0 sample: LR={LR0:.3f}, p={p0:.3f}")

# Under H1
mu_true = 0.7
data1 = laplace_rvs_inverse(rng, size=n, mu=mu_true, b=b_true)
LR1, p1 = laplace_lrt_mu(data1, mu0=mu0)
print(f"H1 sample: LR={LR1:.3f}, p={p1:.3e}")

# Quick type-I check (asymptotic calibration)
m = 1500
pvals = np.empty(m)
for i in range(m):
    d = laplace_rvs_inverse(rng, size=n, mu=mu0, b=b_true)
    _, pvals[i] = laplace_lrt_mu(d, mu0=mu0)
print("approx P(p<0.05) under H0:", float(np.mean(pvals < 0.05)))

H0 sample: LR=0.151, p=0.698
H1 sample: LR=60.434, p=7.608e-15

approx P(p<0.05) under H0: 0.056666666666666664

def soft_threshold(y: np.ndarray, tau: float) -> np.ndarray:
    y = np.asarray(y, dtype=float)
    return np.sign(y) * np.maximum(np.abs(y) - tau, 0.0)


# Bayesian demo: Gaussian likelihood with Laplace prior
# y | theta ~ N(theta, sigma^2),   theta ~ Laplace(0, b)
sigma = 1.0
b_prior = 0.7

# MAP has a closed-form soft-threshold in this 1D model
tau = sigma**2 / b_prior
ys = np.linspace(-4, 4, 401)
theta_map = soft_threshold(ys, tau=tau)

fig = go.Figure()
fig.add_trace(go.Scatter(x=ys, y=theta_map, mode="lines", name="MAP( y )"))
fig.add_trace(go.Scatter(x=ys, y=ys, mode="lines", name="MLE=y", line=dict(dash="dash")))
fig.update_layout(
    title=f"Laplace prior shrinkage (b={b_prior}, sigma={sigma})",
    xaxis_title="observation y",
    yaxis_title="MAP estimate of theta",
    width=850,
    height=420,
)
fig.show()


# Visualize the posterior for a single observation
y0 = 1.2
theta_grid = np.linspace(-5, 5, 2001)

log_like = norm.logpdf(y0, loc=theta_grid, scale=sigma)
log_prior = laplace_logpdf(theta_grid, mu=0.0, b=b_prior)
log_post = log_like + log_prior
log_post -= log_post.max()
post_unnorm = np.exp(log_post)
post = post_unnorm / np.trapz(post_unnorm, theta_grid)

theta_map_grid = float(theta_grid[np.argmax(post)])
theta_map_closed = float(soft_threshold(y0, tau=tau))

fig = go.Figure()
fig.add_trace(go.Scatter(x=theta_grid, y=post, mode="lines", name="posterior"))
fig.add_vline(x=theta_map_grid, line_dash="dash", line_color="black", annotation_text="MAP")
fig.update_layout(
    title=f"Posterior p(theta|y={y0}) with Laplace prior",
    xaxis_title="theta",
    yaxis_title="density",
    width=900,
    height=420,
)
fig.show()

print("MAP (grid)  :", theta_map_grid)
print("MAP (closed):", theta_map_closed)

MAP (grid)  : 0.0
MAP (closed): 0.0

# Generative / differential privacy-style demo: Laplace mechanism for a count query
# For a count, sensitivity Δf = 1. Laplace mechanism adds noise ~ Laplace(0, Δf/ε) = Laplace(0, 1/ε).

true_count = 250
eps_values = [0.25, 0.5, 1.0, 2.0]  # larger ε => less noise
n = 60_000

fig = go.Figure()
for eps in eps_values:
    b_noise = 1.0 / eps
    noise = laplace_rvs_inverse(rng, size=n, mu=0.0, b=b_noise)
    released = true_count + noise
    fig.add_trace(
        go.Histogram(
            x=released,
            histnorm="probability density",
            nbinsx=140,
            opacity=0.35,
            name=f"ε={eps} (b={b_noise:.2f})",
        )
    )

fig.update_layout(
    title="Laplace mechanism: noisy releases of a count",
    xaxis_title="released value",
    yaxis_title="density",
    barmode="overlay",
    width=950,
    height=450,
)
fig.show()


# Tail bound: P(|Noise| > t) = exp(-t/b) = exp(-ε t)
alpha = 0.05
for eps in eps_values:
    b_noise = 1.0 / eps
    t95 = b_noise * np.log(1.0 / alpha)
    print(f"ε={eps:>4}: P(|err|>t)={alpha} at t={t95:.3f}")

ε=0.25: P(|err|>t)=0.05 at t=11.983
ε= 0.5: P(|err|>t)=0.05 at t=5.991
ε= 1.0: P(|err|>t)=0.05 at t=2.996
ε= 2.0: P(|err|>t)=0.05 at t=1.498

11) Pitfalls#

  • Invalid parameters: the scale must satisfy b > 0. Many routines will return NaNs or raise errors if b ≤ 0.

  • Scale vs standard deviation: for Laplace, (\mathrm{sd}=\sqrt{2},b). If you want a target variance (\sigma^2), set (b=\sigma/\sqrt{2}).

  • MGF domain: (M_X(t)) only exists for (|t|<1/b). (The characteristic function exists for all t.)

  • Sampling numerics: inverse-CDF sampling involves (\log U). Clip U away from 0 and 1 to avoid log(0).

  • Degenerate samples: if all observations are identical, the likelihood increases as b → 0+. Any practical MLE needs a small floor.

  • Fitting and medians: for even n, the minimizing set of medians is an interval. np.median returns the midpoint, which is still optimal.

12) Summary#

  • laplace is a continuous distribution on ((-∞,∞)) with location μ and scale b>0.

  • PDF: (f(x)=\frac{1}{2b}\exp(-|x-\mu|/b)); CDF is piecewise exponential.

  • Mean = μ, variance = (2b^2), skewness = 0, excess kurtosis = 3, entropy = (1+\log(2b)).

  • MLE: (\hat\mu) is a sample median and (\hat b) is the mean absolute deviation from (\hat\mu).

  • Sampling (NumPy-only): inverse CDF or difference of exponentials.

  • Common uses: robust modeling (L1 loss), sparse priors (Bayesian Lasso), and privacy noise (Laplace mechanism).

Contents